Section 9 Model Parameterization

These sections are not detailed enough to warrant their own dedicated chapter.

9.1 Parameter Changes

These sections mostly involve changing a few parameter values and are therefore in small sub sections without R-code. Some of this might change later.

9.1.1 Channel parameters revised

This an ongoing unresolved issue #49.

So far nothing has been changed.

9.1.2 Crop parameters verified

The crop database has been updated and stored in the following file:

"model_data/cs10_setup/swat_input_mod/plants.plt"

Discussions to this topic can be found in issue #27.

From the whole database, we only use crops defined in our management files:

Crop management: oats, barl, wwht, swht, pota, fesc

Generic: frst, fesc and others which are minor

The parameters for these crops have been updated to reflect the colder growing conditions.

9.1.3 Soil physical parameters in final form

This has been mentioned in issue #28, and has been closed without discussion. Seems like the parameters are in final form, but no further documentation exists at this point.

9.1.4 Soil chemical parameters in final form

This step will be verified in issue #46. It has been discussed in #19.

The following changes have been implemented, but are subject to change:

Changes made to soilnut1 of nutrients.sol

Parameter	Default	CS10
lab_p	5	20
nitrate	7	8.5
inorgp	3.5	35.1
watersol_p	0.15	0.4

nutrients_sol <- readLines("model_data/cs10_setup/swat_input/nutrients.sol")

header <- nutrients_sol[1]

nutrients_sol_df <- read.table("model_data/cs10_setup/swat_input/nutrients.sol", header = T, skip = 1, sep = "", fill = T, as.is = T, colClasses = "character")

nutrients_sol_df$lab_p <- 20
nutrients_sol_df$nitrate <- 8.5
nutrients_sol_df$inorgp <- 35.1
nutrients_sol_df$watersol_p <- 0.4

write.table(
  nutrients_sol_df,
  "model_data/cs10_setup/swat_input_mod/nutrients.sol",
  quote = F,
  sep = "   ",
  row.names = F, 
)

header_new <- paste(header, "and modified by the CS10 workflow on", Sys.time())

nutrients_sol <- readLines("model_data/cs10_setup/swat_input_mod/nutrients.sol")

writeLines(c(header_new, nutrients_sol), "model_data/cs10_setup/swat_input_mod/nutrients.sol")

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

9.1.5 Impoundment parameters defined

Has been postponed by #20, will be dealt with in #54

9.1.6 Water diversions defined

Deemed as not relevant for this catchment

9.1.7 Point sources parameters added

Is an ongoing issue in #21

9.1.8 Tile drainage parameters defined

Is an ongoing issue in #53.

old_path <- "model_data/cs10_setup/swat_input/tiledrain.str"
new_path <-  "model_data/cs10_setup/swat_input_mod/tiledrain.str"

tiledrain_str <- readLines(old_path)

header <- tiledrain_str[1]

header <- paste(header, "then modified by CS10 workflow at", Sys.time())

tiledrain_df <-
  read.table(
    old_path,
    sep = "",
    header = T,
    skip = 1,
    fill = T,
    colClasses = "character"
  )

tiledrain_df$dp <- 800 # from 1000
tiledrain_df$t_fc <- 12 # from 24
tiledrain_df$lag <- 30 # from 96
tiledrain_df$rad <- 200 # from 100
tiledrain_df$dist <- 8000 # from 30
tiledrain_df$drain <- 40 # from 10
tiledrain_df$pump <- 0 # from 1
tiledrain_df$lat_ksat <- 2 # from 2 (no change)

write.table(header, file = new_path, quote = F, col.names = F, row.names = F)
write.table(
  x = tiledrain_df,
  file = new_path,
  sep = "   ",
  quote = F,
  append = T,
  row.names = F
)

## Warning in write.table(x = tiledrain_df, file = new_path, sep = " ", quote = F,
## : appending column names to file

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

9.1.9 Atmospheric deposition defined

Has been done in section REF

9.1.10 Additional settings verified

Refers to chapter 3.11 in The Protocol and files parameters.bsn codes.bsn Has been discussed in #22

The following changes have been made:

Changes made to codes.bsn

Parameter	Default	CS10
pet	1	2
rte_cha	0	1
cn	0	1
tiledrain	0	1
soil_p	0	1
atmo_dep	a	m

old_path <- "model_data/cs10_setup/swat_input/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"

codes_bsn <- readLines(old_path)

header <- codes_bsn[1]

header <- paste(header, "then modified by CS10 workflow on", Sys.time())

codes_df <- read.table(old_path, skip = 1, header = T, sep = "", 
                       colClasses = "character")

PET method (pet) The OPTAIN project recommends the PET calculation of Pennman-Monteith, which is code 2.

codes_df$pet <- 2 

Channel routing network (rte_cha) Recommended is to start with Muskingum (code 1) and apply variable storage method if it causes problems.

codes_df$rte_cha <- 1

Stream water quality (wq_cha) Recommended to test both for OPTAIN (1/0)

Daily curve number calculation (cn) Code 1 is recommended. (Plant ET)

codes_df$cn <- 1

OPTAIN recommends new version

codes_df$soil_p <- 1

Lapse rate is being tested in REF

Plant growth stress is being testing in REF

Tile drainage

This should be set to 1, but that causes a crash. We will fix this bug in a dedicated secton (REF)

codes_df$tiledrain <- 0

Atmosdep was done in REF and is annual

codes_df$atmo_dep <- "y"

Writing the changes

write.table(header, file = new_path, col.names = F, row.names = F, quote = F) 


write.table(
  codes_df,
  file = new_path,
  col.names = T,
  append = T,
  quote = F,
  sep = "   ",
  row.names = F
)

## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to file

Testing if the setup still works:

simout <- test_swat_mod()
cat(tail(simout), sep = "\n")

9.2 Fixing Tiledrain

Our problem is that some HRUs are drained, but on the wrong soil type. This leads to drains trying to be placed on soils that are not deep enough for them. This was probably caused by bad overlap of the land use and soil maps.

require(mapview)
require(sf)
require(stringr)
require(magrittr)
require(dplyr)
require(purrr)
require(DT)

9.2.1 Finding the problem HRUs

Let us identify where we have this problem. For this we need the soil data, land use map, and landuse.lum file.

basin_shp <- read_sf("model_data/input/shape/cs10_basin.shp")
basin_shp$Id = "CS10 Basin"

lu_lum <- read.table("model_data/cs10_setup/swat_input_mod/landuse.lum", skip = 1, header = T, sep = "")

# need to match IDs 
lu_lum$type <- lu_lum$name %>% str_remove("_drn") %>% str_remove("_lum")

Our drains are set to 80cm depth. Which of our soils are less than that?

NOTE: might need to adjust for drain radius!

usersoil <- readr::read_csv("model_data/input/soil/UserSoil_Krakstad.csv", show_col_types = F)

shallow_soils <- usersoil$SNAM[which(usersoil$SOL_ZMX<800)]

print(shallow_soils)

##  [1] "Mountain"     "CMlen-dy"     "Antropogenic" "Sea_dep_thin" "Humic"       
##  [6] "LVlen-sl"     "End_moraine"  "RGlen-dy-ar"  "RGlep-dy"     "STlen-dy"    
## [11] "STlen-lv-sl"  "STlen-rt-sl"  "STlen-um-sl"

Which of our HRUs use these soils, and are drained?

hru_data <- read.table("model_data/cs10_setup/swat_input/hru-data.hru", skip = 1, header = T, sep = "")

shallow_hrus <- which(hru_data$soil %in% shallow_soils)

drained_hrus <- grepl(x = hru_data$lu_mgt, pattern = "_drn") %>% which()

problem_hrus <- base::intersect(shallow_hrus, drained_hrus)

length(problem_hrus)

## [1] 16

Ok. Just 16 HRUs have an issue. Lets get some info on these 16

problem_hru_data <- hru_data[problem_hrus,]

datatable(problem_hru_data)

hru_shp <- read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
hru_soil <- hru_data %>% select(soil, name)

hru_shp <- left_join(hru_shp, hru_soil, by = "name")

problem_shps <- hru_shp %>% filter(name %in% problem_hru_data$name)

problem_hru_ids <- problem_hru_data$lu_mgt %>% str_remove("_drn_lum")

problem_lus <- hru_shp %>% filter(type %in% problem_hru_ids)


hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_lus$name) == FALSE)
hru_shp <- hru_shp %>% filter((hru_shp$name %in% problem_shps$name) == FALSE)

problem_lus <- problem_lus %>% filter((problem_lus$name %in% problem_shps$name) == FALSE)

Data processing done, lets have a look:

hru_map <-mapview(
  hru_shp,
  map.types = "Esri.WorldImagery",
  color = "green",
  lwd = 1,
  legend = FALSE,
  zcol = "soil",
  alpha.region = 0
)

problem_map <-
  mapview(
    problem_shps,
    color = "red",
    lwd = 3,
    legend = FALSE,
    zcol = "name",
    layer.name = "Problem HRUs",
    alpha.region = 0
  )

problem_lus_map <- 
  mapview(
    problem_lus,
    color = "orange",
    lwd = 1,
    legend = FALSE,
    zcol = "name",
    layer.name = "Problem LUs",
    alpha.region = 0
  )

full_map <- hru_map + problem_lus_map+ problem_map

full_map

Resolution to this problem is being discussed in issue #53

lum_fixed <- read.table("model_data/cs10_setup/swat_input_mod/landuse.lum", skip = 1, sep = "", header = T) %>% tibble()

hru_select <- hru_data %>% select(name, lu_mgt)

hru_select$hru_name <- hru_select$name

hru_select$name <- hru_select$lu_mgt

lum_fixed <- left_join(lum_fixed, hru_select, by = "name")

hru_fixed <- hru_data

9.2.2 Fixing the problem HRUs

hru0218 with landuse a_008f_2 – Remove the drains since it is an isolated LU

lum_fixed$tile[which(lum_fixed$hru_name=="hru0218")] = "null"

hru3637 with landuse a_148f_1 – Remove the drains since it is an isolated LU.

lum_fixed$tile[which(lum_fixed$hru_name=="hru3637")] = "null"

hru3686 of landuse a_115f_1 has the wrong soil type, changing it to “STlv-sl” of its connected field.

hru_fixed$soil[which(hru_fixed$name=="hru3686")] = "STlv-sl"

hru4880 of landuse a_182f_5 has the wrong soil type for the land use. Changing to the “ARdy” of neighboring fields

hru_fixed$soil[which(hru_fixed$name=="hru4880")] = "ARdy"

hru5033 with isolated landuse a_112f_1 can have its drains removed.

lum_fixed$tile[which(lum_fixed$hru_name=="hru5033")] = "null"

hru0692 of isolated land use a_037f_1 can safely has its drains removed

lum_fixed$tile[which(lum_fixed$hru_name=="hru0692")] = "null"

hru0704 has the wrong soil type for its landuse a_037f_2, which should be STlv-sl instead of Sea_dep_thin

hru_fixed$soil[which(hru_fixed$name=="hru0704")] = "STlv-sl"

hru3314 of landuse a_060f_2 has the wrong soil type. It should be STlv-sl, not STlen-rt-sl.

The same is true for hru3376.

hru_fixed$soil[which(hru_fixed$name=="hru3314")] = "STlv-sl"
hru_fixed$soil[which(hru_fixed$name=="hru3376")] = "STlv-sl"

hru2784 (and the within hru3367) is a sizeable field with landuse a_072f_1 and soil RGlen-dy-ar. Tough call here, but best bet is probably to assign soil of the other major field it shares a landuse with (STlv-sl)

hru_fixed$soil[which(hru_fixed$name=="hru2784")] = "STlv-sl"
hru_fixed$soil[which(hru_fixed$name=="hru3367")] = "STlv-sl"

hru3465 of LU a_146f_2 has the wrong soil for the greater field. Changing it from RGlep-dy to STrt-sl

hru_fixed$soil[which(hru_fixed$name=="hru3465")] = "STrt-sl"

hru5214 of isolated land use a_020f_4 can safely have its drains removed.

lum_fixed$tile[which(lum_fixed$hru_name=="hru5214")] = "null"

hru5861 of landuse a_210_f2 does not share the same soil (end moraine) as the rest of its field.

hru_fixed$soil[which(hru_fixed$name=="hru5861")] = "STlv-sl"

hru5926 of landuse a_211f_1 has the wrong soil type (Sea-dep-thin), changing to that of the greater field (STrt-sl)

hru_fixed$soil[which(hru_fixed$name=="hru5926")] = "STlv-sl"

hru2633 of islolated LU a_131_f1 can safely have its drains removed.

lum_fixed$tile[which(lum_fixed$hru_name=="hru2633")] = "null"

9.2.3 Writing modified changes

lum_fixed <- lum_fixed %>% select(-hru_name)
lum_fixed <- lum_fixed %>% select(-lu_mgt)

lum_fixed <- distinct(lum_fixed)

lu_header <- paste("lu table after fixing up by OPTAIN-CS10 workflow on", Sys.Date())

write.table(lu_header, "model_data/cs10_setup/swat_input_mod/landuse.lum", quote = F, row.names = F, col.names = F )

write.table(lum_fixed, file = "model_data/cs10_setup/swat_input_mod/landuse.lum", sep = "\t", quote = F, append = T, row.names = F, col.names = T)

## Warning in write.table(lum_fixed, file =
## "model_data/cs10_setup/swat_input_mod/landuse.lum", : appending column names to
## file

Six drains have been removed.

hru_header <- paste("header header header, delete me and regret it forever! fixed soil types by cs10 worflow on", Sys.Date())

write.table(hru_header, "model_data/cs10_setup/swat_input_mod/hru-data.hru", quote = F, row.names = F, col.names = F )

write.table(hru_fixed, file = "model_data/cs10_setup/swat_input_mod/hru-data.hru", sep = "\t", quote = F, append = T, row.names = F, col.names = T)

## Warning in write.table(hru_fixed, file =
## "model_data/cs10_setup/swat_input_mod/hru-data.hru", : appending column names
## to file

Ten soil types have been changed.

9.2.4 Re-enable tiledrain

Load.

old_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"
new_path <- "model_data/cs10_setup/swat_input_mod/codes.bsn"

codes_bsn <- readLines(old_path)

header <- codes_bsn[1]

header <- paste("modified by CS10 workflow on", Sys.Date())

codes_df <- read.table(old_path, skip = 1, header = T, sep = "", 
                       colClasses = "character")

Re-enable

codes_df$tiledrain <- 1

Write.

write.table(header, file = new_path, col.names = F, row.names = F, quote = F) 

write.table(
  codes_df,
  file = new_path,
  col.names = T,
  append = T,
  quote = F,
  sep = "   ",
  row.names = F
)

## Warning in write.table(codes_df, file = new_path, col.names = T, append = T, :
## appending column names to file

9.2.5 Testing SWAT+

Now we can see if the model runs

source("model_data/code/test_swat_mod.R")

simout <- test_swat_mod()

## [1] "79 of 79 files were copied into the run folder, amounting to a size of  15.77 megabytes"
## [1] "overwriting base files with modified files"
## [1] "8 of 8 modified files were copied into the run folder, amounting to a size of  0.76 megabytes"
## [1] "updating time of simulation"
## [1] "running model for 1 year"

cat(tail(simout), sep = "\n")

##   Original Simulation       12  30  2016 Yr    1 of    1 Time  15:44: 1
##   Original Simulation       12  31  2016 Yr    1 of    1 Time  15:44: 2
## 
##  Execution successfully completed 
## 
##   Date of Sim   6/ 1/2023 Time  15:44: 3

Success.